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ABSTRACT 


Atmospheric turbulence will induce phase and amplitude fluctuations in 
propagating electromagnetic waves, such as a laser beam. Adaptive optical systems 
attempt to compensate for these distortions. The Strehl ratio is a measure of the peak, on- 
axis intensity after propagation through turbulence divided by the peak irradiance for 
vacuum propagation. 

This thesis investigated the probability distribution of the Strehl ratio of a perfect, 
phase-only, adaptive optical system as a function of the atmospheric coherence length, ro, 
divided by the actuator spacing, d. Using an efficient Fourier algorithm and 28 work- 
stations running in parallel, over 850 million computer simulations were performed for 25 
different d/r, ratios in order to produce a histogram of the Strehl ratio probability 
distribution. The results show that the Strehl ratio follows a log-normal probability 
distribution even for very small probabilities. 

A second set of computer simulations introduced intensity scintillation by 
including the log-amplitude variance parameter, c7. Much faster, state-of-the-art 
computer workstations enabled over two billion realizations on 18 machines running in 
parallel for comparable time periods. The trends of these results are more complex and 


will require further research and deeper investigation. 
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I. INTRODUCTION 


A laser beam propagating through the atmosphere will experience degradation due 
to atmospheric turbulence. Temperature differences associated with the turbulence 
generate changes in the air density, and subsequently its index of refraction, with the 
result that parts of the laser beam's spherical wavefront are slowed by differing amounts. 
The field of adaptive optics uses a deformable mirror or a similar device to compensate, 
or correct, for the distortion of the light, effectively reversing the turbulent effects and 
restoring the spherical shape of the wavefront. 

There are three important parameters of significance to this problem. The first is 
the coherence length, or ro, which is a measurable parameter that indicates the severity of 
atmospheric turbulence induced wavefront distortion, and it is extremely variable 
depending on atmospheric conditions at a particular site. The second is the log-amplitude 
variance, o; , which is a measure of the intensity variation of the light perturbed by 
turbulence. 

The third parameter is the Strehl ratio, S, which is a measure of a light beam’s 
intensity at the nominal focus, relative to what it would be if the focusing process were 
diffraction limited. The best one can hope for is to achieve diffraction-limited intensity, 
which can be approached with very good adaptive optics correction of atmospheric 


turbulence effects — but which cannot be reached. 


Pertinent to this thesis is the relationship between r, and o; and the randomly 
varying values of the Strehl ratio of a laser beam degraded by atmospheric turbulence and 
(partially) compensated by adaptive optics. The nominal (or average) value of the Strehl 
ratio will be degraded as ro gets smaller and as o,” gets larger. However the value of the 
Strehl ratio at any instant is a random variable and there is no exact relationship between 
that value and the values of r, and ሮ/. In fact, for given values of r, and or there will be 
many different patterns of turbulence effects on the light beam, each of which will yield 
its own Strehl ratio. A study of the nature of the distribution of the random values of the 
Strehl ratio and the relation of the distribution to the values r, and o,” is the subject of 
interest here. 

The problem of assessing the magnitude and significance of the Strehl ratio 
degradation experimentally is made very difficult by the time and effort required 
collecting these measurements in actual site situations. Not only is the required 
equipment complex and expensive, but also conditions change so rapidly with time and 
one can never “select” particular weather conditions to accumulate data measurements. 
As a consequence a very great deal of time in the field would be required to implement 
such an approach. 

Fortunately, due to exponential advances in computer technology, computer 
simulation modeling is becoming increasingly viable for all fields of endeavor. This 
thesis will simulate the propagation of a laser beam through the atmosphere, with 


selectable values of ro and o, and then evaluate the subsequent Strehl ratio 


measurements. By combining the computer resources across the Naval Postgraduate 
School, I have been able to run as many as twenty-eight Monte Carlo simulations in 
parallel and generate very large sets of the Strehl ratio values from which the probability 
distributions can be calculated with statistically significant results at very low probability 


levels. 
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Il. BACKGROUND 


A. OVERVIEW 


The atmosphere has multiple affects on the propagation of a laser beam. 
Molecular and aerosol scattering and absorption will attenuate the beam. High intensity 
laser beams can heat the atmosphere as it travels — with consequent effects on the 
propagation of the laser beam, which effects are known as thermal blooming. In addition, 
atmospheric turbulence will distort the phase of a beam and make its intensity non- 
uniform as a result of random variation of the atmosphere's index of refraction and multi- 
path interference induced by the turbulence. As this thesis will only address lower power 
laser beams not significantly affected by thermal blooming, atmospheric turbulence will 
be the only factor considered (Greenwood, 1992, pp.3-7). 

Figures (1) and (2) on the following pages show the effect of the atmosphere on 
laser propagation. These figures plot the intensity distribution as a hot-cold differential. 
In Figure (1), the intensity appears as a hot-cold color-scale for a Gaussian beam after | 


long distance propagation through a vacuum. The diameter of the Gaussian waist will be 


proportional in this case to the wavelength divided by the aperture diameter uU » 


However, after passing through a turbulent medium as Figure (2) illustrates, the resulting 


beam pattern includes spreading, intensity and phase variations. 





Figure 1. Laser Intensity Distribution through a Vacuum 


B. COHERENCE LENGTH 

At the heart of describing atmospheric turbulence is the index of refraction 
structure parameter. Denoted by C,, it relates the mean square difference in index of 
refraction at two arbitrary points to the two-thirds power of the distance between the 


points. It is defined by the following equation: (Tatarski, 1961) 


ን (n; -m y) 


Ge u (1) 


273 
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Here, r denotes the distance between the two points and the angled brackets indicate an 


























Figure 2. Laser Intensity Fluctuations through a Turbulent Medium 
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ensemble average value. The z“ dependence arises from Kolmogorov statistics of the 


turbulence. 


In practice, it is usually more convenient to measure C, indirectly, by measuring 


the analogous temperature structure parameter, C} and relating back to C; according to 


the relation 
P 2 


where P is the atmospheric pressure in Pascal and T is the atmospheric temperature in 


Kelvin. Balloon-launched devices with a pair of fast-response temperature probes 


attached are commonly used to measure C7, and from this, C? values are determined. 


Once C? is known along the path of beam propagation, the coherence length, or ro, can 


be calculated through the following equation for a ground to space propagation: (Fried, 


1966) 
0 
2 5 
r, =2.1| 1.46- k? secó jioa | (3) 
0 


In this equation, k=27z//A1s the wavenumber and @1s the zenith angle of the path whose 


length is z. The coherence length is a measure of the electric field auto-correlation length 
and indicates the severity of atmospheric turbulence. It is extremely variable depending 
on atmospheric conditions at a particular site. 

The significance of the coherence length arises with an example from astronomy 
(for light propagating down through the atmosphere), although it applies exactly the same 
for laser beam propagation up through the atmosphere, through the principle of 
reciprocity. When considering the resolution of a telescope, it is a well-understood fact 
that in the absence of turbulence, the angular separation (Qsep) of two stars that can just be 
resolved is proportional to the wavelength divided by the aperture diameter (Born & 
Wolf, 1970, p.415). For conventional telescopes, r, is the diameter of the largest aperture 


that can be used without turbulence significantly degrading the image. The effect is that 


A A 
ወ“ E — —. Asthe turbulence gets stronger, the values for ro are driven smaller and 
Yo 


the resolution limit/separation (sep Will get larger/coarser. For earthbound observatories 
operating at visible wavelengths, ro typically ranges between five and fifteen centimeters 
(Hardy, 1994, p.61). 

In addition to the spreading of an image spot by random phase variations, 
turbulence also induces variations in the intensity pattern. This is a process perhaps most 
simply characterized as “random apodization.” Both phase and intensity variations effect 


the size of the image spot and the system's Strehl ratio. 


€: STREHL RATIO 


There have been literally hundreds of volumes written on the propagation of light 
and its mathematical characteristics, so only a brief summary is provided here, pertinent 
to the problems addressed in this thesis. Figure (3) below represents a coordinate system 


for discussion (Tyson, 1991, p.6). As a means to simplify terms in the exponent of forth- 





Figure 3. Reference Coordinates for Electromagnetic Propagation 


coming expressions, the following normalized coordinates in the focal plane are used: 


_2n(ay _2 (4 
> 5 < 2 (4) and pe E Ol (5) 


For a beam of coherent light at wavelength A and variable amplitude function, 4e 4^? , 
where € =7+ig can be complex, the intensity of light at a point P on the image plane a 
distance z away is given by: (Born & Wolf, 1970, p.461) 

| 


The notation used here is illustrated in Figure (3), where the variables (p,0,z) constitute 


2 


(6) 


Aa? 
AR? 





KP) = | “ሚወ 
0 0 








a cylindrical coordinate system. R 1s the slant range from the center of the pupil to point 
P, o represents the deviation in-phase from a perfect sphere about the origin of the focal 
plane, and 77 represents the intensity deviations. 

From the standpoint of adaptive optics, and pertinent to this thesis, the most 
important quantity in the above expression is the ¢ term, commonly referred to as the 
wavefront phase perturbation. If no aberrations are present (diffraction-limited case), then 
the intensity is a maximum on-axis (r=0), also called the Gaussian image point: (Born & 


Wolf, 1970, pp.461-2) 





2 
Aa’ 
I io 5 5 2 nte | (1) 


The Strehl ratio, S, is a measure of a beam's intensity at focus (or in the far-field) 


relative to what it would be if there were no optical system induced or turbulence induced 
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aberrations. It is defined as the ratio between the intensity at focus or on-axis in the far- 
field of an aberrated beam and the intensity of an unaberrated beam at the same position. 


The remaining aberrations are represented by ¢, and the Strehl ratio can be expressed as: 


(Born & Wolf, 1970, pp.462) 


I( P. 27 2 
S(v = 0) em) „I | “ህሠሠ 
0 0 


2 
I i-o 25 


: (8) 








For phase only adaptive optics, intensity variations are ignored (because they 
cannot be corrected), and 77 1s often set to zero. If both the intensity and phase effects of 
the beam at the pupil are unaberrated (or dp = zero), the Strehl ratio will reduce to unity. 
S = 1, which is to say that the intensity at the focal plane is diffraction-limited. This is 
the unattainable goal of perfect adaptive optics; to correct the wavefront in such a manner 


that dp = 0 and the Strehl ratio 1s unity. 


D. COMPUTER SIMULATION OF BEAM PROPAGATION 


A word needs to be said concerning the implementation of beam propagation 
using computer simulation. There are essentially two methods. The first method is most 
intuitive to the casual observer, by defining a beam profile and actually propagating it by 
using Fourier transform matrices as the Optical Transfer Functions (OTF's). One can 
create the random phase screen (representing the turbulence), then take the Fourier 
transform and multiply by the input beam pattern. The propagation distance 15 important, 


because the underlying principle is Huygens-Fresnel propagation. This method of 


simulation can be extremely accurate, although multiple passes are frequently necessary 
to represent long propagation lengths through different turbulent areas (Davis, 1994, pp. 
34-36). One other important point is the fact that both amplitude and phase effects from 
the turbulence occur with each propagation step. For the second method, the amplitude 
effects must be characterized separately from the phase effects, and then combined at a 
later time by use of convolution (Fried, Sep 1997, p.ii). 

The second method is based on reciprocity in optical propagation. It allows us, 
instead of considering propagation of a laser beam to the far-field through multiple layers 
of turbulence — instead to consider a distorted laser beam arriving at an aperture after 
propagation through turbulence and then to simply consider the propagation of that laser 
beam, in one step, to focus. Use of this second method allows us to make use of our 
knowledge of the statistics of the beam's distortion after propagation through the 
distributed turbulence to avoid the multi-step calculation of propagation through 
distributed turbulence. Ignoring intensity variations, which we do in this section (though 
not in the next), we simulate the phase variations of the laser beam after it has propagated 
through the distributed turbulence by a single phase screen. (To have simulated 
propagation through the distributed turbulence, we would have needed a number of such 
— weaker — phase screens.) 

The simulation of the phase screen requires knowledge of the power spectral 
density of the medium through which the beam will propagate. This context requires an 
accurate power spectral density for the phase screens. Fortunately, much research has 


been completed in this area by pioneers in the field, such as Tatarski, Kolmogorov, and 
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Fred. Without showing the full development, it can be shown that since the phase 
structure function follows a five-thirds power law, then the power spectral density, ®, , 


has the form: (Fried, Sep 1997, p.2) 


ወ,(፳)= ,”” መ:”” . (9) 
where 
5/6 -] 
A- (2) 1B = 0.04579117422. (10) 


This equation for the power spectral density provides a means for generating phase 
screens that contain the proper atmospheric phase statistics. This formulation allows easy 
scalability to different turbulence conditions through the r,-term. 

Section III of this thesis includes the effects of the turbulence along the 
propagation path in terms of the wavefront distortion at the optics aperture, without the 
effects of intensity scintillation included. The value of r, characterizes the wavefront 
distortion and the phase structure has a five-thirds power-law dependence (Fried, 
Sep.1997, p.1). The efficiency of this method is greater than the first method described 
above, because the entire propagation can be represented by one phase screen, requiring 
only one Fourier transform to complete the simulated propagation. Additionally by 
convolution, a two-dimensional phase screen (256 x 256) contains separate realizations. 

The effects of intensity scintillation can be simulated separately and in the same 
manner as the phase effects, which is to say finding an expression to represent the power 


spectral density for the intensity scintillation pattern. Intensity scintillation is an effect 


that cannot be corrected by any conventional form of adaptive optics (Fried, May 1997, p. 
li). This would be extremely troublesome except for the fact that intensity scintillation 
effects in the atmosphere are much smaller than the phase effects (Greenwood, 1979, p. 
550). Hence, adaptive optic systems are proving very effective by just compensating for 
the phase distortions. Nevertheless, the effects of intensity scintillation are non-negligible 


and will be addressed, later in Chapter IV. 


Ill. STREHL DEGRADATION DUE TO PHASE EFFECTS 


A. OVERVIEW 

These investigations utilized computer simulation source code written by Dr. 
David Fried in the MATLAB programming language, and as such, required modifications 
before applying to this particular problem. In the course of accomplishing the 
modifications, a considerable amount of effort was required to understand the intricate 
programming techniques used to arrive at an efficient algorithm. After all, the ultimate 
task is to produce computer simulation code that can produce very large samples of 
random strehl ratio measurements efficiently in order to investigate very small probability 


events. 


An example here is appropriate: Consider a turbulence strength of interest being 
an average value (rp = 10 cm), and one measurement of the strehl ratio from an adaptive 
optics compensated system results to S = 0.35. Is this average? Is it extremely good and 
luck was on our side? Or is it much worse than had been hoped, and luck definitely was 
not on our side? One very good way to answer this question 1s to collect a large number 
of sample measurements of the strehl ratio when the rp value is always 10 cm, and then 
look at the average and standard deviation of the sample. The larger the number of 
samples, the more confident one can be in saying that, for example, the 0.35 strehl ratio 


was indeed an average and there is a significant probability that the next measurement 


will be as poor as 0.25, ... or not. It is the efficient development of very large numbers 


of sample values that is the core of the work reported here. 


B. PHASE SCREEN GENERATION 


To generate a random phase screen, one must first generate a suitable ‘random’ 
starting set of numbers, and then scale to the desired statistics. Remembering that the 
spatial frequency components are continuous, a substitution will be made to account for 
discrete representation in the computer simulation environment. Discrete spatial 
frequencies will be noted by ‘m °’ in place of the continuous ‘x ’ spatial frequency, both 
in units of cycles/unit length. As such, ‘ w(m di will be defined to be the N x N array of 
random values — each value with its corresponding one of the N x N spatial frequencies. 
In the phase screen generation algorithm that 1s used here, the values of w(m ) are chosen 
from a Gaussian distribution with zero mean and unity variance — in other words, 
standard white noise. The values of w(m) will be complex: (Fried, 1996, p.3) 

w(m)= w,(m) + iw, (m)., (11) 
with w, (m) and w, (m) each being real Gaussian random variables with unity variance. 

As stated previously, from equations (9) and (10), the power spectral density 
produces a suitable random phase screen. That is to say that the phase screens generated 
will have a five-thirds power law structure function. For this to be so, from equation (9), 


the spatial frequency components of the power spectral density must be proportional to 


the minus eleven-thirds power of the magnitude of the spatial frequency (Fried, 1996, 
pp.2-3). Accordingly, we write 


-11/3 


¢,(m)«,| |m| (12) 


where, because the power spectral density will be proportional to the square of the 
magnitude, to arrive with a minus eleven-thirds power law, the square root has been 
applied: (Fried, 1996, p.4) 
هين‎ elm Em). me a3 
0, f |m|-0 

Here, the a-coefficient encompasses the ability to scale the ‘strength’ of the phase screen 
so that it will correspond to the desired value ro, and will be discussed in the next 
paragraph. Also, realizing that the zero spatial frequency component will only affect the 
average (or dc) value of the simulated phase screen it has been set to zero, avoiding an 
infinite result. 

The value of the a-coefficient is developed and given by Fried to be dependent on 
the size, N, of the array: (Fried, 1996, pp.13-15) 


ie ENS 


Nu Ner (14) 
53 (N) - CN) 


a= / 


where u = 0.98965794 , and f — 6.8838772, and the values for S2(N) and S)(N) are given 
in Table (1). 
The random array # (m ) conforming to the power spectral density characteristics 


of a random turbulent phase perturbation, can be converted back into the spatial domain 


I 


2.27148846e-07 2.26348792e-05 
3.49154895e-09 5.239791 46e-07 
5.43330535e-11 1.15360543e-08 


8.48088179e-13 2.46158319e-10 
1.32479989e-14 5.14359846e-12 
2.06986786e-16 1.05908706e-13 
3.23411700e-18 2.15760798e-15 





Table 1. Screen size (N) dependent a-coefficient scaling factors. 


. ላ * ው == =l =- . . . 
via the inverse Fourier transform: ø, (ic) > م‎ (F), in the continuous case — and, in 


the analogous discrete case: 9, (m)— > ¢, (p). This will be accomplished in the 


discrete two-dimensional case, according to the following equation: (Fried, Sep 97, p. 4) 
2ጫ” ኀ p 
Pp(p)=N" > me rim E). (15) 
m,n=] 


where m =(m,n). Once in this spatial domain form, @, (p) will consist of a complex 


random array, the real and imaginary parts of which each will correspond to a random 
phase pattern. This phase pattern is suitable for high spatial frequencies but has a 
deficiency in the lowest spatial frequency components. Dr. Fried states, 

This is manifested by a structure function for this pattern which has the 

appropriate five-thirds power law dependence for separations much 

smaller than the real-world physical size of the screen, L, but deviates 

very significantly from that power law behavior for separations as large as 

L/10, and even decreases in value as the separation increases beyond L/2 

(Fried, Sep 97, pp.4-5). 


In order to compensate for this deficiency in the lower spatial frequency 


components, Fried has developed a suitably chosen random tilt term that can be added to 
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Pp (p). The results with this added term prove sufficient, and the resultant simulated 


phase screen will be accurate, i.e. its structure function will follow a five-thirds power 
law for all separations. The extra tilt-term can be seen developed and validated in the 
reference source (Fried, 1996). However, in this treatise, this additional tilt term is not 
necessary since the adaptive optics will have no problem canceling that part of the phase 
variation. This particular exercise will be limited to laser apertures under the correction 
of adaptive optics, and as such, the adaptive optics will introduce their own transfer 
function — eliminating most the wavefront distortion, to be discussed later. The pertinent 
point here, is that the adaptive optics can correct the tilt terms in the distorted wavefront 
very effectively, thereby eliminating the need to include them in the first place (Fried, 


1998). 


6 ADAPTIVE OPTICS COMPENSATION 

As mentioned previously, the adaptive optics component will introduce its own 
Transfer Function (TF), and as such, will require treatment in this simulation. To be 
precise, the adaptive optics will correct all phase fluctuations to the extent that the size of 
the sub-aperture will allow. For example, if it were possible to build an adaptive optics 
system with infinitely small sub-apertures, then the wavefront distortion could be 
corrected perfectly. Allowing for finite size sub-apertures there will be residual errors 
after adaptive optics correction; we shall consider adaptive optics systems with finite size 
sub-apertures, thus introducing the deviation from perfection whose statistics are studied 


here. In addition, we will ignore latency effects caused by the time delay required by the 
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adaptive optical system. The imperfection of the compensation provided by the adaptive 
optics is represented by what 1s called the residual phase error, by the deviation of the 
corrected wavefront from perfection; this is called fitting error. 

In this simulation, the residual fitting error will be developed from the phase 
screen through the use of its transfer function mentioned earlier. As developed in 
Greenwood, 1979, the adaptive optics corrected residual phase error can be represented 


by the following equation; returning to a familiar continuous treatment, we write 


Pro (7) = Jexp(2x ik -7) T ¿0 (K) bp (i) dk , (16a) 
d ا‎ 9 aol") ar (16b) 
= 10 (K)Qp(K) 


where T, (X) is the adaptive optics fitting error transfer function written in frequency 


terms. $ p(x) is the Fourier transfer of ¢,(7), which is the phase without correction from 
adaptive optics, as previously discussed. As shown in Greenwood, 1979, the value of 


T o () is given by the following equation: 
2J (a| z |d 4J, (z| K |d 


-| a| K d -| 7% ld |. (17) 


where d denotes the effective size of the sub-aperture components in the adaptive optics. 





1 40 (<) 5 





As can be seen from Figure (4), the 7,,(X) transfer function has the effect of attenuating 


the lower frequency terms and passing the high frequency terms, analogous to a high-pass 


filter in electronics (Greenwood, 1979). 
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Figure 4. Greenwood Transfer Function Parallels a High-Pass Filter 


D. STREHL COMPUTATION 

The computation of the Strehl ratio turns out to be the most time-intensive portion 
of the simulation. However, through the use of a fast Fourier transform (FFT) and clever 
assumptions, one statistically random phase screen realization can generate many 
simulation results. In the case of one phase screen represented by equation (15) and 
simulating the random perturbations a light ray has encountered through the turbulent 


medium, the Strehl ratio can be computed by comparing the intensity calculation with the 
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phase screen to that intensity calculation of a reference intensity with zero phase 
perturbations (diffraction-limited). 

For an aperture whose center is at 7 , the signal power, P(7), comes directly from 
squaring the signal strength, S (r). The signal strength calculation is most conveniently 


expressed in the following manner, which facilitates future development: 





S(F) = jal? = | exp(id, (P) dr”, (18) 


where ¢,(7') represents the previously discussed random phase distortion function 


evaluated at r'. The aperture function, A(x), is defined as follows: 


1, if |x|<+ 


(19) 
0, if otherwise 


10- 
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Thus 4 | represents an aperture of diameter D on the r' — plane, with its center at 


r =r. The representation provided by equation (18) is beneficial two-fold: a) the 
reference intensity is easily calculated up-front by replacing the ø,(F') with ¢, [ሆ') = 
zero, and b) the calculation of the signal strength with phase distortion and adaptive 
optics compensation can be completed easily by replacing the ¢$, (7 | with 
办 人 ( ( = ¢,0(F') from equation (16). 

It is here that Fourier transforms can be used to provide a tremendous increase in 
productivity. As described in the previous section and also as represented in equation 


(18), lt is assumed that the 7 vector is positioned at the center of the array (although not 


explicitly stated as such) and that the 7’ vector moves with the integration along all of the 


22 


points over the aperture. In fact, the 7 vector does not need be positioned at the center, 
and as Fourier transforms will allow, the 7 vector can be positioned at every point and 
continue with the integrations through the virtual use of the two-dimensional convolution. 


Following Fried, 1997, the phase exponential will be defined: 
U(r) = exp ፉሬ ሆ)ኦ (20) 


so that equation (18) becomes: 
r-r 
S(F) = |4| —— | U(r’) dr". 2] 
(F) Al > | (7) (21) 


Written this way it is manifest that S(7) is equal to the convolution of A with U. A 
convolution indicated by the usual (*) designation, of the two elements in the integrand of 
equation (21) would appear as follows: 


— f — 
KE oA 





se-4 UE), (22) 


where the S (r ) result would now be an N x N matrix of the same dimensions as the U (F ') 


matrix. Following the treatment shown by Fried, and independently applying convolution 


manipulation, equation (22) is transformed again using Fourier transforms: 

S(F) = [exp(2zik-7)A (<) 0) dE, (23) 
where A (<)and U (x) are defined by the following equations: 

A(R) = | exp(-22ik-7) A(F) dF. (24a) 


Ü (r) = | exp(- 221K -F) U(r) dF. (24b) 


m 
U) 


One other point of note concerns the value of the Fourier transformed elements, 


and in particular the value of A(&). As opposed to the case of U (Z), which includes 


calculating a two-dimensional FFT of the randomly generated U (F ), the value of 4 (X) 15 


a closed-form result depending on the geometry of the aperture. In this case, a uniform 


circular aperture 1s used, and the corresponding result can be expressed as follows: 


(25) 


dg) - on D xem 


4 3 |D 
This expression allows a more refined representation in the computer code, avoiding the 
computer rounding inaccuracies that arise during the calculation of the discrete FFT. A 


smooth illustration of the Fourier transform of a circular aperture is shown in Figure (5). 


E. INCORPORATION INTO COMPUTER CODE 

All of the previous sections contain the mathematical background for efficiently 
generating many simulations of Strehl ratio calculations for a given set of initial 
conditions. As will be detailed later, the different exponential degrees to which ‘many’ 
can cover will depend on the speed of the computer being used as well as how many 
computers (or processors) can effectively be utilized in parallel. This section will present 
the applicable formulations to take into account discrete implementation as opposed to 


the continuous treatment previously discussed. 
The last discrete formulation presented was Equation (15) representing the 28 (p) 


phase perturbation term in frequency space. The next two equations present the discrete 


24 


1600 
1400 


1200 





1000 
800 
600 


400 





Figure 5. Fourier Transform of a Circular Aperture (Spatial Frequency Domain) 


version of the adaptive optics compensation, T,,, (x ) from Equation (16): 


x rd, a IN and (26) 


(27) 





Incorporating Equation (27) into the resultant phase perturbation and performing the 


inverse discrete Fourier transform result in the following two equations: 


$us lm)- T, (2) >, 95) exp 201%, E), and (28 


人 人 | ፈር 2) (29) 


m,n=1 
The culminating incorporation into the discrete adaptive optics compensating phase 


Screen M (p) , 1s given by the anti-climactic expression: 


U ,o(P)=explig,o(P)} (30) 
Its discrete Fourier transformation follows: 
ም N p 
U o (m) = > al) cxp| 2%, Z) (31) 
p,g=1 N 
In order to discretely represent the Fourier transform of the aperture, a conversion 


of Equation (25) yields the following two discrete representations: 


X =rD\R+R, and (32) 


A(m) “ወ ES (33) 


Possessing discrete formulations for both U,,(m) and A (m), the discrete conversion of 


Equation (23) follows easily: 


S(p)- N^? y Ü (m) A(m) zi 2zik,, 2) (34) 


m,n=] 

After completion of Equation (34) in the computer code, we are effectively left 
with a statistically random N x N matrix of signal strength values, all for the same 
aperture diameter D and the same r, value. The signal strength values can then be 
squared and divided by the reference intensity value, yielding an N x N matrix of Strehl 


ratio values. Figure (6) is a block diagram representation of the functional flow 
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throughout the computer code. It starts from the very beginning of variable initialization 
and reference intensity computation through the finishing — which is only truly finished 
when the desired number of Monte Carlo realizations has been reached. Numbers, for 
instance through the completion of this thesis, range from two loops in the process of 
testing code accuracy to 83,000 loops in the preparation of ten billion realizations for one 

particular r, value. A later version of the computer code, written in Matlab and named 
‘Sin.m’*, is provided in Appendix (A). As such, this version includes extra code to 


account for intensity scintillation effects to be discussed later in chapter IV. 


F. RESULTS 

1. Preparation 

Over the two-week span of 12-26 March 1998, results were gathered from 
implementation of the aforementioned computer code, which at the time, was named 
‘turro.m.’ The following parameters will be analyzed as a result of those computer runs, 
although all may be changed to run any particularly interesting case. An aperture diameter 
size of 1.2 meters will represent the telescope, with adaptive optics sub-aperture 
components having an effective diameter size of 15 centimeters. The resultant D/d ratio 
in this case becomes 8.0. The simulated size of the phase perturbation screen will be 3.0 
meters, thereby fully covering the simulated aperture. Finally, the range of eleven r, 


values will be applied, from 12.6 to 14.1 centimeters in increments of 0.015 centimeters. 


¢ ‘Sin’ is short for ‘scintillation’ in this particular case, not to be confused with any religious meaning. 
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Figure 6. Block Diagram of Functional Flow in the Computer Code 


The resultant d/r, ratios in this case become 1.190 down to 1.064. These parameters are 


summarized in Table (2). 
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Aperture Diameter | 2000) 
Sub-Aperture Diameter | 0.15 (m) 
D/d Ratio 8.0 
Screen Diameter 3.0 (m) 
Values ofr, (meters) 0.1410, 0.1395, 0.1380, 0.1365, 0.1350, 0.1335, 0.1320, 

0.1305, 0.1290, 0.1275, 0.1260 

1.0638, 1.0753, 1.0870, 1.0989, 1.1111, 1.1236, 1.1364, 
1.1494, 1.1628, 1.1765, 1.1905 

Values of Scintillation Zero 
Screen Size(N) 250 X256 


Table 2. Initialized Parameters for Results of 'Turro' Computer Runs, 12-26 March 1998. 



























Resultant d/r, Ratios 











As mentioned earlier, the design of the computer code allows many simulations to 
be run on different computers, effectively parallel processing without actually linking the 
processors. In fact, the storing of the values in a sparse matrix (histogram fashion), 
allows the programmer to gather the results from each separate computer that was utilized 
and add all of the results together into one combined sparse matrix. Figure (7) below 
displays the results that were generated from the above initialized values. The 
observation of “830 million per line” indicates the number of Monte Carlo simulation 
results that built each apparently Gaussian shaped bell curve. Each run through the main 
loop of the program yields 2 x 256° or 131,072 realization results. One quick calculation 
reveals that 6336 runs through the main loop were needed for each r, value, yielding a 
total of 69696 runs for all eleven values of turbulence severity. 

To appreciate the magnitude of this number, consider the following actual 
numbers obtained during the generation of the lines on Figure (7). One run through the 


main loop of the computer code, executed in Matlab on a 200 MHz Intel® Pentium® Pro 
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Figure 7. 12-26 March 1998 ‘Turro’ Computer Code Results 
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(circa 1997) personal computer (PC), requires approximately one minute. The required 
run time would result to 484 days which is clearly well over a year for completion. 
Utilizing access to a multitude of Sun® Sparc®-workstations (circa 1990) across the 
Naval Postgraduate School campus proved to be the saving grace. Although the time to 
complete the main loop on any individual Sun was significantly slower than the PC, the 
greater number was invaluable. 

A familiarity with the Unix® operating system is essential to efficiently gather 
data from multiple runs across the campus network. One needs first to check ‘the 


prospective machine to determine the processes already running and then start the Matlab 


programs running in the background before moving to the next machine. The “rsh” 
command is needed to log-on quickly, complete business on the machine, and then exit to 
“rsh” on another. At the time of execution in March of 1998, there were seventeen 
Sparc®-stations in Ingersol Hall, thirteen in Root Hall, and eight in Spanagel Hall for a 
total of 38 prospective machines. Only 28 of the machines turned out to be ‘available’ 
after checking to ensure that nobody else was running their own programs, taking up 
precious processor usage time. Over the course of nearly a week on 28 machines, the 830 


million simulations were realized. 


2: Examination of Results 

The data displayed in Figure (7) was gathered in such a way as to allow 
examination in the following manner. In order to analyze low-probability events, the data 
must be converted from a histogram showing a bell curve into a corresponding 
probability distribution graph. The histogram is essentially a probability density curve for 
each curve. Computing the cumulative probability is necessary in order to generate the 
next graph, displayed as Figure (8). A curious circumstance of working in such extreme 
low probabilities is the number of ‘nines’ (high), or ‘zeros’ (low). A probability of 
0.0000001 = 107 is displayed in the graph’s y-axis as (-7), short for 1.0E-7, and 
conversely a probability of 0.9999999 = 1-10” is displayed as (7), short for 1.0E-7 or 7- 
nines (high). 

The next order of business is that of prediction --- or at the very least, 1s there 


some natural mathematical law that is discernable from the data? Referring to Figure (8) 
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Figure 8. Strehl vs. Probability with Crossing Line 


again, there is a vertical broken line drawn at Strehl = 0.30. One can visualize an infinite 
number of parallel curves above and below the eleven curves shown, and the vertical line 
graphically cuts through all eleven and extends outward into the infinite number of other 
such curves not shown, nor even generated. Changing the axes of the graph and plotting 
differently, Figure (9) below shows five such vertical lines plotted as d/r, versus 
probability. The lines appear to be very straight. Nevertheless, exactly how straight are 
they? 

To check the straightness of the lines, it is necessary to generate a separate set of 


data at less severe turbulence values, or higher r, values. Take the higher Strehl = 0.370 
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Figure 9. Turbulence Severity vs. Probability for selected constant Strehl lines 


line for example from Figure (9). One can “fit” an exact line to the data points that are 
plotted as an apparent straight line. That same Strehl = 0.370 fitted line can then be 
extrapolated further to the bottom-left of the graph. Figure (10) displays that detailed 
extrapolation for the top three lines from the previous Figure (9). Figure (10) also 
displays the results, shown as small circles according to the legend, of Strehl probability 
data generated separately according to the initialized variables shown in Table (3). The 
initialization values shown in Table (3) were chosen to provide “checker” data points. As 
plainly seen from Figure (10), the “check” data, labeled as “test” data, does not fall on the 


extrapolated line from the “fit-to” data. At the time, this was a very discouraging result. 
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Aperture Diameter 1.2 (m) 
Sub-Aperture Diameter | 0.15 (m) 
D/d Ratio 8.0 
Screen Diameter 3.0 (m) 


0.1590, 0.1560, 0.1530, 0.1500, 0.1470, 0.1440 
0.9434, 0.9615, 0.9804, 1.0000, 1.0204, 1.0417 


Values of ro (meters) 
Resultant d/r, Ratios 
Values of Scintillation | Zero compensation 
Screen Size(N) 256 x 256 


Table 3. Initialized Parameters for “Test” Data of 'Turro' Computer Runs, 12-26 March 1998. 






Evaluation that is more extensive led to the following question. If not straight on a linear 


scale, how straight is the line on a log scale? 
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Figure 10. Strehl Extrapolation Lines on a Linear Scale 
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Keeping the dependent axis the same, but changing the independent axis to a log 
scale proved to be the answer. Figure (11) shows the results from such a change. The 
fitted line was fitted to the logarithm, base 10, of the "fit-to" data. The same 
extrapolation was performed and in the same manner, the small circles represent where 
the "checker" data resulted. As can also be seen plainly from Figure (11), the data fall on 


a very straight line, indeed, albeit a straight logarithmic line. 
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Figure 11. Strehl Extrapolation Lines on a Logarithmic Scale 
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IV. STREHL DEGRADATION DUE TO PHASE AND 
SCINTILLATION EFFECTS 


A. OVERVIEW 

It is a fact that adaptive optic systems are indeed proving very effective by just 
compensating for phase effect distortions along the path of light propagation. However, it 
is also considered interesting and possibly extremely helpful to model the effects of 
intensity scintillation and combine the model with that already discussed in an attempt to 
evaluate the same low probability Strehl events. As mentioned previously, the effects of 
intensity scintillation can be simulated separately and in the same manner as the phase 
effects, that of finding an expression to represent the power spectral density. Once that 
expression 1s determined and then validated, it can be included into the ‘turro.m’ program 
to again generate many simulations for a given rọ value and a given value of intensity 
scintillation, o,’. 

For the following development, it is assumed that the particular optical system in 
question has the capability to perfectly correct phase-only perturbations with its adaptive 
optics components. We will continue to use the principle of reciprocity to analyze a laser 
beam propagating outward and upward though the atmosphere towards a target in space. 
And as before, we will only be concerned with an evaluation of the statistics of this 
scintillation effect. As such, and following Fried (from Fried, 1997, p.3), the statistical 


notation, (F), called the turbulence induced random log-amplitude at position r, will 


represent all of the effects of turbulence that will be of concern to us in this thesis. The 


next section will address the details of its development. 


B. SCINTILLATION SCREEN INCORPORATION 
The purpose of this section is to develop a similar random screen to the phase 
perturbation screen from the previous chapter, with this particular screen, d,, having the 
statistical characteristics of atmospheric turbulence induced intensity scintillation. More 
explicitly, we wish to able to express an optical field, U (F ), in the following expression: 
U (7) - 42 exp(9, (7)). G5) 
This will prove desirable as we can again use equation (21), developed earlier in chapter 


three and reproduced here for the benefit of the reader: 





S(F) = jaf U(F) dr’. en 


1. Log-Amplitude Fluctuation Statistics 


The random log-amplitude, /(r), as developed in the turbulence literature from 
Kolmogorov to the present, may be considered to be a random function characterized by a 
mean value / and a covariance function C, (7). These two functions are defined by the 
following equations: (Fried, 1997, p.5) 

)36( ريحت 


ad C(F)=([e(F+17)-2][e(7-4)-2)) . (37) 
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A very instrumental parameter follows, called the log-amplitude variance and 
denoted by c. It corresponds to the 7 =0 value of C, (7), which is to write explicitly: 
TIC (38) 
The importance of this parameter stems from the following relationship that can be 
shown from conservation of energy principles: 
Deo, (39) 
This states simply that given a value of c^, we empirically know the mean value of the 
random log-amplitude function (Fried, 1997, p.5). 

The following equation, (40), follows from more detailed development presented 
in Fried, 1997. The pertinent points to be taken from this development include the fact 
that the C,(r)function depends on the path and the particular wavelength. The 
propagation path in this thesis runs from s=0 on the ground to s => œ at the space end 


of the path. The notation k again refers to the optical wavenumber related to the 


wavelength by the expression k=27/4. Presented without derivation, we write: (Fried, 
1997, p.6) 


co = Bk’ lig) se ds, (40) 


PATH 


where the value of B is given as B x 0.56315761. Hence, we are left with a means to 


; : : 2 
measure, or in this case, generate a given value of o. 


The other important parameter, the covariance function C, (7), is needed in order 


to generate the power spectral density, b. Once given a power spectral density, we can 


impart that statistic onto the random generated values representing the spatial frequencies 
— exactly as done in the previous chapter when producing the phase perturbation screens. 
The important results will again be presented, allowing the interested reader to review the 
source document for details. After appreciable development in Fried (1997, pp.6-7), we 
write: 

C,(F) = lexp(2z is -7) ወ,(፳) dF, (41) 
where it is noted that the two functions are Fourier transform pairs. The power spectral 


density, ®,, is given by the equation: 


$, (c) - AF Qa) P || | ct i-es (ante) | ds, (42) 


PATH 


where the value of A is given as A = 0.65150030 and the path is the same as above. 


2 Screen Formulation 

Incorporation into scintillation screen computer code follows very closely to that 
presented in the previous chapter. The generation of the array of random values is 
accomplished in the same manner, by defining the spatial frequency components of the 
array, w(m). These elements will again be complex and take on random values with zero 
mean and unity variance. For the sake of clarity, the two components of the m vectors 


will be x,, and «, , with the values before randomization initialized as follows: 
= M 
K, and, s 7 : (43) 


where m will be defined: 
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e if l«m«iN (44) 


Lon-N-D if iN«m«N 

Ihe second step will also be the same as the random values are scaled to be 
proportional to the magnitude of the square root of the log-amplitude's power spectral 
density. Corresponding to equation (12) in the earlier development, this can be expressed 


in the spatial frequency domain as follows: 


im) = win) (a.[®,m)) , (45) 


where the a scaling factor is defined as œ = N?/L in order to account for the physical 


size of the screen. This will be apparent in the next step, where the result of the inverse 
Fourier transform, as realized by an inverse FFT, will now produce a suitable magnitude 
result. (Fried, 1997, pp.10-11) 

The third step corresponds to the recently developed equation (40), replacing the 


integration with a summation: 


(p)2N^ y L(m) expl 22%. (46) 


m,n=1 
For the sake of clarity, the 2(p) result is a two-dimensional complex array of random 
values denoting spatial coordinates, remembering prior to the final step that the mean 
value here is still zero. 
All that is left to do in the fourth and final step is to apply the mean value as 
defined in equations (36) and (39). This is realized by computing equation (40) in 
discrete increments (i.e. estimating the integration by selecting a finite number of 


intervals in the path). The mean value adjustment is made, yielding the screen(s) and 


4] 


written in the following two relationships: 
60D), 2 81j-o; , (47) 
and $ DB), =3 Up) j 0. (48) 

It is at this point we have a statistically sound representation of the log-amplitude 
function, and the interested readers are referred to the source document (Fried, 1997, pp. 
11-15) for verification of the accuracy with this method. These screens represented by 
equations (47) and (48) are now ready to implement in the computer code as an extension 
to equation (30), written as follows: 

U or: (p)= exp Li dao (p) j-exp 2 (p)j (42) 
There is no deviation in the remainder of the ‘turro’ computer from this point on. 

A word about the real-world practical significance of the turbulence induced log- 
amplitude value is relevant here. The unit of measurement for o, is the square Napier, 
which is denoted as Np*. Actual turbulent airflow measurements range from an ideal of 
zero through a typical measurement of o,’ = 0.05 Np’, up to a severe turbulent condition 
of o» = 0.15 Np’. Jt is with this o, value that the turbulence severity will be represented, 


just as the rọ value represented the turbulence severity in the phase perturbation of the 


beam. 


C. RESULTS 
Over the two-week span from 21 October- 4 November 1998, results were 


gathered from implementation of the combined computer code, which at this time 1s 
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named ‘sin.m.’ The following parameters will be analyzed as a result of these computer 
runs, although just as earlier mentioned, all initialized parameters may be changed to run 
any particularly interesting case. An aperture diameter size of 1.1 meters will represent 
the telescope, with adaptive optics sub-aperture components having an effective diameter 
size of 10 centimeters. The resultant D/d ratio in this case becomes 11.0. The simulated 
size of the phase perturbation screen will be 3.0 meters, which again will fully cover the 
simulated aperture. 

With additional code to account for the intensity scintillation, the average time to 
complete one main loop of the code is increased to approximately one minute thirty 
seconds on the same 200 MHz PC. Consequently, the range of r, values will be reduced 
to three, namely 10.2, 11.1, and 12.0 centimeters. The resultant d/r, ratios in this case 
become 0.980 down to 0.833. Finally, each of the r, values will be coupled with five 
intensity scintillation values ranging from a very severe o, = 0.11 Np? all the way down 


to zero scintillation. These parameters are summarized as well below in Table (4). 


Aperture Diameter 1.1 (m) 
Sub-Aperture Diameter 0.10 (m) 
D/d Ratio 11.0 
Screen Diameter 3.0 (m) 


Wavelength (A) 0.840 (microns) 

Values of ro (in meters) 0.1200, 0.1110, 0.1020 
Resultant d/r, Ratios 0.8333, 0.9009, 0.9804 
Values of Scintillation (in Np’) 0.11, 0.07, 0.03, 0.01, 0.00 

Screen Size(N) 256 x 256 





Table 4. Initialized Parameters for Results of 'Sin' Computer Runs, 20-31 October 1998. 
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The long time delay between the two sets of runs (seven months) proved to be 
beneficial in one respect: it allowed new workstation computers to be acquired, greatly 
improving processor speeds and memory availability. The following histogram results, 
displayed in Figure (12) below, represent two billion simulated realizations of the 
combined computer code. The results are also displayed with a legend, clearly illustrating 
the gradual Strehl ratio decline for each of the three r;-values as the intensity scintillation 
gets increasingly worse. 

The long time delay between the two runs also proved to be detrimental in another 
respect: the results were actually gathered after the author's graduation date and thereby 


represented a successful stopping point in the life of the thesis. Time is not available in 
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Figure 12. 20-31 October 1998 ‘Sin’ Computer Code Results 
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this thesis for detailed examination of the results, as we were able to do in the previous 
chapter. Only cursory remarks will be made during the following presentation of the 
various 'snapshot figures, with the intention of stimulating thought for further 
examination in another thesis, possibly. The cursory examination will be presented in 
two segments. The first segment is shown as “blow-up” views of each of the three 
separate ro-value groups, followed by the five views of each of the o,”-value groups. 

Close evaluation of Figure (12) reveals one very evident feature: the bell-shaped 
characteristic of each curve gets wider for each increase in 0,. An even clearer 
manifestation of this is seen in Figures (13), (14) and (15) where it can be seen that the 
slope of the curves decreases as o,” increases. The rotational differences in the slope of 
the probability lines towards the clockwise direction are indicative of the widening. 
There is also an appearance that the lines would actually cross high and to the right, if in 
fact, the lines were extended further to the right. 

We will make two remarks here concerning the previous observation. The first 
remark is a confirmation of earlier hypothesis and very limited actual measurement. The 
G -value has a greater impact on the variance of the Strehl ratio measurement. The 
decrease in slope as the intensity scintillation increases, so indicated in the three previous 
figures, is in direct support of the increase in variance. Statistical calculations on each of 
the bell curves support an increase in both the standard deviation and the variance as the 


intensity scintillation increases. 
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Figure 13. Strehl vs. Probability Lines at d/r, = 0.833 


The second remark will require extra work in the future for validation. This 
concerns the overall r,-value after both r, and o,’ parameters have been combined. For 
example: we first selected an r,-value (i.e., ru. = 12 cm) and then selected a scintillation 
value (c^ = 0.05 Np5), without readjusting the r,-value for the worsening that will occur 
as a result of the addition 5 scintillation. In fact, the actual ro-value would worsen (to 
approximately ry = 10.5 cm) because the added intensity scintillation has an effect on the 
rg- value measured at the aperture. An approximate solution would be to centroid align a 
set of intensity variance curves leaving d/r, constant. Unfortunately, since the shapes of 
the histograms change slightly for different d/r,, this is not exact. Further research will be 


needed to resolve this interaction. 
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Figure 15. Strehl vs. Probability Lines at d/r, = 0.980 
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It is the author's contention that as the intensity scintillation increases, the starting 
ro- value would need to be weaker (1.e., r, = 13.5 cm), so that after applying the same 
scintillation value (o, = 0.05 Np’), the desired ro-value would deteriorate (to ru. = 12 cm). 
This ‘fix’ would transpose each of the lines to the right in proportion to the increase in 
c^, forcing the intersection to move down to the center position of 0.5 probability — 
effectively having all five of the lines rotate around one central position. 

The remaining five figures (16—20) present the same information organized in a 


different manner. In examining these figures, it will be noted that the Strehl ratio shifts to 


the left as each figure represents a constant o,’-value with varying phase-only ro-values. 
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Figure 16. Strehl vs. Probability Lines for cM = 0.11 Np! 


48 


2.06 x 10? per curve 


D/d=11.0 1 7 


o? = 0.07 Np? / / 
f A 
/ / 


Probability 


0.4 0.45 


Strehl 


2.06 x 10? per curve 


ከከከ 
ሀዩ = 0.03 Np? 


Probability 


0.4 05 
Strehl 





Figure 18. Strehl vs. Probability Lines for ©,” = 0.03 Np’ 
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Figure 20. Strehl vs. Probability Lines for o, = 0.0 Np? 
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V. CONCLUDING REMARKS 


A. LESSONS LEARNED 
The list of the many lessons that were learned will be shortened to include only 
the most important ones. The big lessons that must be learned very well include the 


following: 


Become intimately familiar with the Matlab programming language. 


O Learn the basics of the Unix OS, well enough to move files around (using ftp 
and telnet), and execute files in the background. 


(] Get access to the fastest workstation processors available that can run 
Matlab. At the time of this printing, on the NPS grounds, this includes the SGI 
Octane® and the Sun Ultra 60 Sparcstation®. 

Smaller lessons are important as well, because they will often be the cause of 
pounding headaches. First, ıt is ımperative to have one different program for each 
separate machine, in order to save to separate data files. The separate data files will 
continue to grow as the program progresses, and when finished, the data files can be 
loaded one-by-one into Matlab and recombined by adding together the individual sparse 
matrices. Second, the random number generator must be well understood. A kev 
statement that should be inserted into each file is the following: 
“randn(‘state',sum(100*clock)).” Otherwise, Matlab will begin each individual session 


with the exact same ‘state’ and then follow with the same random numbers every time -- 


it is known as a pseudo-random number generator. 


Sil 


During the course of this thesis, several problems were encountered while 
attempting to run Matlab scripts in the background on different workstations around the 
NPS grounds. First of all, the Matlab code must be a script to run in the background. 
Functions will not process. Good information can be found on this topic at the Matlab 
web site. Next, it 1s highly recommended to check the processes running on a particular 
machine — before starting a Matlab process in the background. The ‘ps’ command will 
perform this function, with the ‘ps -aux’ switches working well on SunOS systems, and 
‘ps -aef * working well on Solaris and IRIX systems. Anything less than 99 percent 
processor usage time 1s detrimental to the speed of processing one's program, as well as 
increasing the chances of death by interference. Working over long weekends and holiday 


breaks is a good practice to follow. 


B. POSSIBLE FUTURE RESEARCH 

This would best be stated as ^where do we go from here?" There are several 
avenues of research that can be pursued. Just as the scintillation was added to the 
simulation from the initial effort, more effects can be modeled and effectively added to 
the simulation. Optimization of the computer code so that it may be run on faster 
mainframe computers is an interesting field of work as well. This section will address 
just a few of these avenues, and rest assured, close liaison with Dr. Fried is strongly 
advised (vital) for any work surrounding this computer code. This is especially true now 


that it has been changed quite a bit by the author. 
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Within the operational field of adaptive optics, there are two especially promising 
topics thought to be “doable” with respect to realistic simulation on computers. The first 
topic 1s called branch points. Branch points are manifested in actual field measurements 
from adaptive optics systems, specifically, when attempting to correct phase perturbations 
in the presence of strong scintillation. Very simply, adaptive optic corrections are 
‘measured’ by sensing the phase of adjacent light rays in the focal plane, and then 
calculating differences between the measurements. When there is zero intensity in one 
particular spot on the focal plane, there is a 27 anomaly in the phase which the adaptive 
optics system cannot cope with properly. This is known as a branch point from 
mathematics theory and it causes a malfunction in the sensor. Dr. Fried has some 
thoughts on how this might be represented in this type of computer simulation. 

One other promising topic is called /atency within the adaptive optics field. 
Inherent in the process of phase corrections is the fact that a correction is inevitably being 
applied based on a measurement taken a short time (on the order of milliseconds) prior. It 
should be possible to develop an algorithm capable of making a more informed 
correction, based not only previous data, but also an educated (algorithmic) prediction. 
Successful efforts in this area would reduce feedback problems associated with the 
sensors and enable higher Strehl ratio capabilities. 

There are three more topics that present challenging difficulties, but not much 
thought has been devoted to solutions. They are listed briefly here, intended only to 


stimulate thought: 


(] Measurement noise, 
O Non linear signal processing, and 


O Servo effects not covered by a time delay. 


With respect to computer code manipulation, the need for speed in processing 
time never diminishes. Short of waiting for the next generation faster processor and 
faster memory, the only other way 1s to "streamline" the code. During the course of this 
thesis, we attempted to use a Matlab compiler with the intention of generating an 
executable file in C or C++. Conversion into executable code would enable adaptation 
onto much faster mainframe computers, such as SGI Cray machines, for example. 
Unfortunately, there were two important functions within Matlab (‘sparse’ and ‘spline’) 
not supported by the Matlab compiler — curious as that was. We also considered 
converting the entire program into C++, but quickly found that daunting task more suited 
to a team of C++ programming specialists. Recently, newer versions of the Matlab 
compiler are being released with claims of support for previously unsupported functions, 
and would thereby be worth trying. 

One other important point concerning Matlab that may be important when porting 
to multi-processor machines is the limitation that it cannot utilize them efficiently. 
Fortunately for this application, the parallel nature of the code not only allows multiple 
runs to process on separate machines, but multiple runs capable on the same multi- 


processor machine. Each of the processors within the computer will take dedicated time 
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for its respective program, and there is no loss in processing time — effectively getting 


twice the output in the same amount of time on a dual-processor machine. 


C. CONCLUSION 

In this thesis, over two billion computer-simulated realizations verified a log- 
normal relationship between the Strehl ratio and its probability distribution value for a 
given d/r, value, and similarly a log-normal relationship between the d/r, value and its 
probability distribution value for a given Strehl ratio value. Extremely low probabilities, 
down to 1.0 x 10°, correspond to the number of realizations collected, with Stability in 
the lines down to a probability of 5.5 x 10°. Most recently, over the course of a week, 21 
computer workstations plus another 5 dual processor workstations (effectively acting as 
10) were running 31 separate copies of the computer code to generate 15 curves. 

This computer code is available to run simulations for any given adaptive optics 
aperture arrangement in any given atmospheric turbulence condition. Through the 
principle of reciprocity, the code is applicable to laser beam transmission up through the 
atmosphere, as well as light collection in an astronomical telescope. It can be scaled to 
run any desired number of realizations, given the proper computer resources and realistic 
time constraints. The histogram results are conducive to probability analysis and can 
therefore be used in conjunction with other criteria to make informed decisions about 


prospective adaptive optics system design 








APPENDIX A. MATLAB COMPUTER CODE ‘SIN.M’ 


Sin.m 

TU AAA AAA de 4 40 2 7 | 131 1 11111111222 T 

ረ 

4 Scripts called Sin (short for scintillation), to incorporate Dr. Fried's 


/ programs. Both of these function programs are designed to 

/ run using variables referenced to an r-sub-D (ro) which is equal to the 
% $12e1 La of the phase screen. i-e-1 D is actually incorporated as D/ro. 
0 


Z This script vill utilize the same calculations and fix the aperture 
4 Sizes. D4 to a fixed sub-aperture sizes d+ and run through a range of 
4^ (ro) values. The 'fixed' parameters, Di di and Lı therefore must 

% be scaled for each run through different (ro) values- 


% The scintillation has been incorporated also. to run through a range of 
4 values using sigma-sub-L squared. 


Z 

AIA I R R E de Tode I Ihe da وك لد ولد‎ E ያ ዳያ ጸም O Ia e መሪ ጠዳ መ IL Te Te Ie Ta ee Ie Ie de ILI IE 
7 

% l1) N = Lb thru 1024 possible choices (in powers of 2) 


። ሮን) Required Matlab functions or scripts for this to work at compile time: 
7 -- Scinfun.m 


7 -- avint.m 
ረ -- hv57.m 

^ -- psdgen.m 
Y 


% 3) Choose desired values for: 

7 -> ro [range of values]; 

% -- sigma [range of values]; 

y -- D. aperture diameter in meters: 

A -- d4 sub-apertures effective diameter. in metersi 
7 -- L^ effective length of phase screen, in metersi 
Z -- دمع‎ number of loops for desired computation 

-- filename, for saving data 


LLL IL II eI Lo d e, ILL መ ም ም ውቅ ቁ ውብ ም ም ም ም ም ቁማ ም ዳዳሪ TT 
tO = clocks 
% Define parameters for Ro manipulation 


N = 25b; 


ro = £.040,..037:.034J; 
sigma = E.11..07..03..0125 


= 1.15 % aperture diameter. in meters 
= 0-15 % sub-aperture diameter, in meters 
= 35 „ screen length. in meters 


lengthíro)5 


D 
d 
L 
aa 
bb = lengthísigma)5 


% Scale the physical parameters to the varying Ro value 


EES-OE./roi 
DD = D./ro5 
dd = d-/ros 


clear L D di 


E 


% Calculated results! for S-functions in denominator of alpha 

S1=f2.2?14664be-07%, 3.u815u885e-DH4 5-43330535e-11. 8.480881l79e-13 
L-32479969e-14. 2-0b98b ?6be-lb. 3-23411700e-14835 

Se=Ce-2b348?Fee-05. §-239791L4be-O7. 1-153b0543e-08. 2-4b1583154e-10 
5 . د12-عطلا8 5 35 لاز‎ 1-0590870be-13ı 2.157b0798e-151I5 

sl=Sl(log2(N)=3) 3 

ፍሮፎ=ኣሮ(109ሮ(ከክ)-3) 35 

clear Sl S25 


M = N/23 
alpha = sqrtí(Lb.0883877?x((2xM)aA (-5/3)-D.7987558x(2xM)A (-2))/(s2-s1))5 


/ Define the spatial frequency range for phase screen generation later 
f=EO: (M-1) (=M): (=1) 13 

fD=f.AD: 

Fe=ones(2xM,1)xf05 

Fe=Fe+Fe' 4 

Fe(1.1)=eps5 

clear f fe; 


F3 = alpha * (F2.A(-11/12))5 
F3(1,1)=05 
clear alpha 


ppuUtzrl1 
IntenRef-E135 


for qq = l:aa 


% Generate the spatial frequency scaled by 1/LL and then in 2-2 
kappa = TO: ((N/2)-1) (-N/2): (-1)1/LL(qg)5 

kappa = meshgrid(kappa) 34 

kappamag = sqrtíkappa-Ad+kappa' -AD)5 

kappa = linspace(min(min(kappamag) ) 1max (max (kappamag) )12*N)4 


^ Generate the spatial freq. transform of the apertures which again has 
4 been scaled by the varying Ro value 

x = pi*kappaxDD(qq)3 

We = D.25xpix(DD(qq)Aa2c)x(L[1 2xbesselj(l1=x(2:(2xN)))./x(2:(2xN))D)5 


% Calculate the reference intensities for Strehl computation later 
ppuUtD = spline(kappa.Wt) 3 
U = sqrt(e)xones(N)4 % Phase screena U 
Ut የየሂፎ(ሀ)5 % Phase screen transform, Ut 
人 
። 


Ut ppval(ppUtüskappamag?5 Aperture transform, Wt 
V = abs(ifftecUt.*Ut))4 Uniform reference intensity for diff-limit 


/Concatenate results for later multiple use 
IntenRef = CIntenRef. V(Clal)aeds 
ppWt = Cpplwts ppwtdO J; 


end 
clear x Wt ppwWtO U Ut ل‎ 5 


Prob=sparse(1,1.,0.,aa*bb;,10e3,10e3*aaxbb) 3 
uu = 1:1e45 


randn('state' .sum(1D0xclock))5 
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^ Begin major looping through ranges of (ro) and (sigma) values 
for fn = 1:750; 


count = 1; 4 Keep track through the interior loops. 
for nn = 1:asıh Loop through the (ro) values 
for mm = 1:bb5 % Loop through the (sigma) values 


scrn-fn = ዖሦልበበከ(ሮችቨ-4ሮ#ቹቨ)+፥3 ችዮ”3ኋበገ።በበ(ሮቹቨ ፡,ሮቶሽቨ) 5 
sern=fn = F3.xscrn=fn5 


x=pi*saqrt(F2) *dd(nn)/LL (nn) 4 
x(1.1)=eps5 
T=sqrt(1-(2xbesselj(1l,x)+./x).a2-(Yxbesseljl21x)./x).A2)5 


clear x3 
T(1,1)=03 
scrn = scrn-fn.*T3 clear T5 


scrn = iffto(scrn) * ((LL(nn))Aa(5/b))5 


Function call to get the scintillation screen scaled to sigma value 
[scrnScJI = Scinfun(Nisigma(mm) ) 4 


kappa = ED: ((N/2)-1) (-N/2): (-1)1/LL (nn)>5 
kappa = meshgridíkappa)5 
kappamag = sgrtíkappa-Ad+kappa' -4d)5 


U = sart(e)*lexplit*treal(scrn)).*exp((real(scrnSc)-sigma(mm))))5 
ffte(U); 

ppval (ppWt(nni:).kappamag) 4 

V - abs(ifftg(Ut-.xUt))3 


= 
ct 
ዘ wu 


prob = round(10De3*(V(:).Aad/(IntenRef(nn)))): 
ii = findíprob == 0)5 

prob(ii) = 15 

Prob = Prob + sparse(countiprobiliaa*bb,10e3) 54 


U = sart(2)*lexp(i*imag(scrn)).*exp((real(scrnSc)-sigma(mm))))4 
Ut = ffte(U): 
V = abs(C(ifftga(Ut.*Ut))3 


prob = round(10e3x(V(:).A2/(IntenRefínn))))5 
ii = find(prob == 0)3 

prob(ii) = 13 

Prob = Prob + sparse(count,probil.raa*bb.l0e3) 


count = countths 


end 
end 


% Save the data every few runs in case of unsolicited stopping. 
if mod(fni20)==0 
ttt=etime(clock,tQ) 34 
save Sin?Dl Prob fn ttti 
)م015‎ [ ‘frame #' intestr(fn) '.' J ); 
end 
end 
clear nn scrn-fn scrn Fe F3 ppwt kappamagi 
clear U Ut Ut V kappa prob scrnSci 


D9 


fn = €'done ' intestr( fn ' 33 


ttt=etime(clock,tQ)3 
uu uu/leus 


4^ Save the data here at the end 
save Sin?01l Prob fn ttt uui 


% Execute TimeDisplay M-file. if it exists. 
if exist('TimeDisplay.m') 

TimeDisplay(ttt)5 
end 


4 Important to have this 'quit' when running on 
4 Silicon Graphics SGI machines 
quit 
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APPENDIX B. REQUIRED MATLAB FUNCTIONS FOR ‘SIN.M’ 


Scinfun.m 
2111# * ة[ [ | | |[ |[ |[ |[ [|[|[ |[ | [|[ |[|[|[|[|[|[|[|[ »| [ | | |[ |[ | [| |[ |[ |[ | [ |[ |[ | |[ |[ [ |[ |[ ز ز<زذ<‎ 000017 0 0 


Function: called Scinfun,. to incorporate Dr. Fried's two programs, 
psdgen and log-amp-screen- As the adaptive optics are phase-correcting 
only. no corrections will be introduced here. The turbulence will be 
varied through the severity of the C-sub-N squared. which in turn will 
be set to the values used in the Turro script runs- This will also 
require a wavelength selection, not needed in Turro. 


NNNNNNNNWN 


7 ببببببدد----+-+-+ 2 - 112 11 | * 11 
7 

% N = 16 thru 1024 possible choices (in powers of ( 

17 

VAAL Ala VELA VU Ee IL UL Ve IEEE VI VL eae eU, መመሪ ታረ መመረ ታሪ... 


function Escrnasiglsql -» Scinfun(Nsro) 


4 Define parameters for Ro manipulation 


lamb = 0.8Ue-b:i / 840 nanometers - near infra-red light propagation 

D = 0-45 4 physical aperture diameter (m) 

L = 35 4^ physical phase screen diameter (m) 

4 Scale the physical parameters to the varying Ro value 

LL = L/ros / all arrays (1 x length(ro)) 

DD = D/ros 

clear Di % keep (ro) for later conversion to actual value (cm) 


% Generate the PSD for the scintillation in the turbulence 
air = logspace(0O,logl10(3e4) 1250) ° 
CN8D » hv5?(air)35 


Z Still need to coordinate the ro with CNe -- basically dial the 
% CNe value using the ro value as the knob 


% We will assume that (ro) from hv5? equals 5.7 
hveeS? = -05?5 


randn('state’,sum(100*clock))5 


((roxL)/hvee5?)a(-5/3)5 
aaa * CNeDs35 


aaa 
CNe 


] 5 1 9 15 0 ضع ول 5 م>21 26512 مم د‎ ) 31-1 ነ 


% Define the spatial frequency range for phase screen 
mt = 1:(N/2)5 

mt Emt-1 (((N/2)+1):N>)-N-115 

kappa = mt/LL5 

kappa = meshgrid(kappa) 4 

kappamag = sqrt(kappa-ac + kappa’ -Ac) 4 

clear kappa? 


ዘ 
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srPSD = saqrt(ppval(ppPSD.ikappamag) )3 
clear kappamagı 


alfa » N^ec/LL3 

scrn-fn = randn(N)+1i*randn(N)4 
scrn-fn = (alfaxsrPSD).x*scrn-fn5 
clear srPSD3 


sern = iffte(sern-fn)5 


avint.m 


LALALA AL LLL AL ISAM ily ل‎ aT 


NNN NNNNNNN 


Function, called avint(x.iy.ixlo.ixup) to perform integration 

of a curve defined by data points (x.y) using overlapping 

parabolas. Integrates from xlo to xup. These do not need to be one of 
the actual abscissa points. x(i)- Must have xlo < xup. These default 
to xlo = x(1) and xup = x(length(x))5 


Coded from Davis 2 Rabinowitzw "Methods of Numerical Integration," 
pg-483 => t3b 8/1/91 


function I=avint(xay.ixlo.4xup) 


n= 


length(x)5 


if nargin<3., xlo=x(l)4xup=x(n)5i end 


syl=xloi 


for ib=]:n 


if x(ib)<=xlo, break-ı end 


end 
lb=minímax(e=ib).-n-21)5 


for j=n:-1:1 


if xup>=x(j). break. end 


end 


j= 


max(ibamin(jain-1)-1)3 


sum=D5 
for jm=ib:j 


x1=x(jm-1)31x2=x(jm)5x3=x(jm+1)5 
tl=y(jm-1)/((x1l1-x2)x(x1-x3))15 
1ር=ህ(1)/( (ደር LIA 2 MD JA 

ት A AE A O E 
A=tl+td+t35 

B=-(x2+x3)*tl - (x1l+x3)xt0 - (xl1+x0)xt35 
C-x8xx3*tl * xlxx3xteo + x1xxcxt31 


if jm==1b 

ca=A;cb=Bicc=C34 
else 

ca=.5*(A+ca)icb=-5*(B+cb) 5cc=-5*(C+cCc) 5 
end 


sum=sum + cax(xea3d=syla3)/3 + cbx.5x(x2a2c-sylac) * ccx(xe-syl)3 
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ca=A3icb=Bicc=C5 
syl=xd5 
end 


I=sum + ca*r(xupa3-syla3)/3 + cbx*.5x(xupae=syla2) + coxtxup-syl); 





ንበር m 

Veo IS TIU TIR Tee de 1e Ve TE de 10 0070 2L 109. VII, d LIL LL ht dd AI ገረዳ.፡፡፡ብ፡፡፡፡፡፡፡ ፡፡፡፡ 2222272222 SA 
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% Function, called hv5?, computes the HVS/? Cnae profile - This was 

% converted to m from the Fortran source in the t0SC library. 

1 

^ alt is an array of values and the function returns an array of values. 
7 

7 MODEL = XXXO? 

7 

Z modelname = 'Modified Hufnagel-Valley’ 

7 

TOM 7772 LILIA Ai ASA ANNAN 


function cn2 = hv5?(alt) 
model = 2107; 
cne » Dxalti 
inatmos = (0 <= alt <= 24000): 


if all(~inatmos) 
return 
end 


al = alt(inatmos); 

v= ceil(model/100); 

ene(linatmos) = (2.de-53xal-AlD .x ((v/2?)a2xexp(-a1/1000)) 
* l.e-lb*exp(-al1/15000)2x28.?7?* 

low = (O <= alt < 11000); 


ርበፎ(10ህ) = cn2e(low)+1.%e-14xexp(-alt(low)/100)5 


psdgen.m 
2.፡ገ፡ገ፡፡፡፡6556.::፡፡፡/፡፡:29 መመመ መረረ ሚማ የ የረ ረረረሥረሪገ2፡.. .. ፡: . .. 


Function, called psdgen. used to generate results for the PSD of the 
log-amplitude covariance; Phi-sub-l(k). and for the log-amplitude 
variances. sigma-sub-l-squared. It takes as inputs values for the wave 
length. and the value of C-sub-N-squared at various positions: s; 
along the propagation path- The PSD result is reported out in spline 
form, suitable for use by the function ppval. 


oe on TN DNT EN N as Ne 


--دبببببب7ببببببب7ب7بببببنببب-بببب-بب ب 2030 


INPUTS 
S Clim) = Position along propagation path (in meters). 
CN? (lam) - Value of C-sub-N-squared at positons along the 


propagation path (in m4(-28/3)). 
Uave length (in meters). 


NNN NNN 


lambda (141) 


% 

% OUTPUTS 

% siglsq (1.1) = Log-amplitude variances sigma-sub-l-squared 

7 (in Npad). 

7 ppPSD (145m) - Spline coefficients for the evaluation of the 

% power spectral density---at a spatial frequency 

% whose magnitude is given in cycles/meter--- 

7 Cin Npae/(cycles/m)ac). 

7 

LAMA A e Vs ILL و7 يك‎ IIa TM LA n e IIT AL A ሪፈሪራታታሪሪታረሪ ረረሪሪሪ .፡ ረሪሪ tae 


( 5 12-1 ) + 5 ) مع و0 5 م لل 9 285مم +50 5191 1 function‏ 


k = 2xpi/lambdas 
8) ዳዓ :> 
B = 2a(-4/3)xpia(3/2)/(3xcos(pi/l2)*gamma(2/3))5 


siglsq = B*ka(?/b)* feval(FisiCNe-*s-a(S/b)) 5 


A = (2a(5/3)xsqrt(pi)/39)x(gamma(11/b)/gamma(10/b))5 
kappamax = 20*(k/(e*xpixmax(s)))5 

kappa = kappamax*logspace (-5:0:333)5 

kappaml13 = Axka2x(2xpi)a(-8/3)xkappa-A(-11/3)5 
kappasq = (2*pi)ackkappa-ae/ki 


PSD = zeros(1:333)3 
for n=1:333 

PSD(n)= kappam1l13(n)*feval(FisiCNe-.*(l-cos (kappasq(n) *s)))5 
end 


ppPSD = spline(kappa.PSD) 5 
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APPENDIX C. EXTRA INTERESTING RESULTS 


Along the course of completing this thesis, not everything went according to plan, 
nor as one would want to happen. The following mistake was the cause of nearly a 
month in lost time, but hopefully at a cost of deeper understanding. At some point during 
the preparation phase of the combined phase and scintillation program, “sin.m,” the 
author inadvertently input unrealistic values for the aperture diameter and its sub-aperture 
diameter. The “unnoticed until examination of results” mistake led to an unrealistic D/d 
ratio of 4.0 and unintentionally broadened the Strehl range of the Gaussian-shaped bell 
curves. The resultant "straight" cumulative probability lines were noticeably not straight, 
and, in fact, the “sin.m” program could not duplicate the previous “turro.m” results when 
zero scintillation was initialized. After a line-by-line check of the two programs, the 
different aperture initializations were noticed, but not before one billion simulation 
results were developed and examined. The results from that data set, called at the time 
“combine.m,” are reproduced here in the appendix for the perusal of interested readers. 

As presented in the body of the thesis, the following Table C-1 1s a concise 
summary of the initialized values for this set of results. It can be seen that we ran five 
values of r, versus seven total values of E (including the value of zero scintillation). 

The following Figures (C1-C3) should be self-explanatory after the thesis 
presentation for the subsequent runs of “Sin,” although we will mention two points of 


interest. First, the range of Strehl ratio realizations is quite large, from approximately 
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0.25 up to 0.75. This can be physically explained as having effectively only four sub- 


apertures 


Aperture Diameter 0.40 (m) 
Sub-Aperture Diameter 0.10 (m) 
D/d Ratio 4.0 
Screen Diameter 3.0 (m) 


Wavelength (A) 0.840 (microns) 
Values of ry (in meters) 0.1260, 0.1200, 0.1140, 0.1080, 0.1020 
Resultant d/r, Ratios 0.7937, 0.8333, 0.8772, 0.9259, 0.9804 
Values of Scintillation (in Np’) 0.10, 0.08, 0.06, 0.03, 0.02, 0.015, 0.00 
Screen Size(N) 256 x 256 





Table C-1. Initialized Parameters for Results of 'Combine' Computer Runs, 10-26 September 1998. 


across the entire aperture, not able to provide much correction at all. Second, as a 
consequence of such a wide spread in Strehl ratio realizations, the probability lines in 
Figure (C2) are very long and not as straight as noticed in all of the other runs (which 
were done with a more realistic D/d ratio in every case). This can be physically explained 


as approaching the linear boundaries of zero and one on either side of the Strehl axis. 
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Strehl Distribution 


Strehl Bin Width = 0.0001 
107x10. per curve 3 
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Figure C2. Strehl vs. Probability Lines at d/r, = 0.794 
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1.07 x 10? per curve 
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Figure C3. Strehl vs. Probability Lines for 6,” = 0.06 Np” 
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